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Abstract 

We propose a two fluid theory to model a dilute polymer solution 
assuming that it consists of two phases, polymer and solvent, with two 
distinct macroscopic velocities. The solvent phase velocity is governed 
by the macroscopic Navier-Stokes equations with the addition of a force 
term describing the interaction between the two phases. The polymer 
phase is described on the mesoscopic level using a dumbbell model and its 
macroscopic velocity is obtained through averaging. We start by writing 
down the full phase-space distribution function for the dumbbells and 
then obtain the inertialess limits for the Fokker-Planck equation and for 
the averaged friction force acting between the phases from a rigorous 
asymptotic analysis. The resulting equations are relevant to the modelling 
of strongly non-homogeneous flows, while the standard kinetic model is 
recovered in the locally homogeneous case. 



1 Introduction 

The importance and continued research interest in the non-homogeneous flows 
of suspensions and polymeric fluids may be attributed to the richness of the fluid 
mechanical phenomena often associated with such flows (such as migration, slip 
and the Fahraeus-Lindqvist effect |10| ) as well as their prevalence in polymer 
processing, oil recovery, bior neology and microfluidic devices, for example. In 
this paper we shall be concerned with the rigorous mathematical derivation of 
the governing equations for non-homogeneous flows of dilute solutions of FENE- 
type dumbbells. 

There is some agreement on the macroscopic equations describing the rhe- 
ology of non-homogeneous flows of dilute polymer solutions. The issue of the 
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compatibility of three different formulations of these equations was addressed 
by Beris and Mavrantzas [4] in 1994. The governing equations, specifically those 
for the elastic stress tensor and the polymer number density, were obtained in 
the references cited by the authors using an non- homogeneous kinetic theory [5] , 
a continuum two-fluid Hamiltonian model |18| and a body tensor formulation in 
continuum mechanics [53] • Beris and Mavrantzas showed that a stress diffusion 
term that appeared in the constitutive equations derived by Bhave et al. [5] and 
Ottinger also arose in their two-fluid Hamiltonian model, provided that a 
particular dissipation term was kept in their derivation and not omitted as in 
the original paper |18j . A term that was quadratic in the stress appeared in 
the constitutive equation of Ottinger [23J but not in the other two formulations. 
This term, however, would not lead to qualitative changes in the predictions of 
the three models, at least for small values of the velocity gradient. We note here 
that the derivation using a kinetic theory model of a closed-form constitutive 
equation for Hookean dumbbells by Bhave et al. [5] was only possible because 
of their use of the homogeneous form of the Kramers expression (Eqn. (40) of 
their paper). After redeveloping and correcting the analysis of Bhave et al. [5], 
Beris and Mavrantzas showed that the polymer number density equation was 
the same for all three formulations. 

An important phenomenon associated with nonhomogeneous flows of poly- 
mer solutions is that of molecular migration. By this we mean that the centre 
of mass of a molecule does not move at all times with the same velocity as that 
of the surrounding solvent at the position of the centre of mass. The argument 
of Tirrell and Malone [29], among others, that a macromolecule will try to min- 
imize its free energy by moving to a region where stresses are lower was rejected 
by Aubert and Tirell [2] and Aubert et al. [3J who concluded that there was no 
mechanism for coupling the free energy driving force to the motion of the centre 
of mass. Aubert and Tirrell [5] showed that for a free-draining dilute solution of 
bead-spring chains where the solvent velocity v s could be expressed as a second 
order Maclaurin series 

v s (x,t) = v(0,t) + x- V« s (0,i) + ^xx : VVv(0,i), 

the only migration of molecules could be in directions where v s has non-zero 
Cartesian components. In particular, cross-streamline migration was stated to 
never occur in rectilinear flows (i.e. parallel flows with only one Cartesian com- 
ponent) but could occur in curvilinear flows. The authors then illustrated their 
point by showing that in Poiseuille flow of a dilute polymer solution no mi- 
gration occurred in the radial direction, although the polymer lagged behind 
the solvent along streamlines. In circular Couette flow, however, it was demon- 
strated that cross-stream migration took place. The analysis of Aubert and 
Tirrell [5] was generalised, with similar conclusions, by Sekhon et al. [57] and 
Brunn and Chi [5]. These authors showed, however, that the introduction of 
hydrodynamic interaction made cross-stream migration possible even in parallel 
flows. A similar point was made more recently by Ma and Graham |17| who 
stated that in the absence of hydrodynamic interactions between polymer seg- 
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ments no migration is found in shear flows without streamline curvature. It 
should be noted here, however, that wall effects on the cross-stream gradient 
of polymer stress were not considered in arriving at these conclusions. In the 
work of Beris and Mavrantzas [3] , Bhave et al. [5] and Ottinger [53] , discussed 
above, cross-stream migration was possible in wall-bounded non-homogeneous 
flows because of a term proportional to the divergence of the polymer conforma- 
tions in the definition of the polymer mass flux, the conformation tensor being 
anisotropic due to the presence of the walls. In the Fokker-Planck-based ap- 
proach of Lozinski et al. [TB] to the problem of plane Poiseuille flow of a dilute 
polymer solution, similar considerations of wall effects on the allowable poly- 
mer conformations led to a non-homogeneous flow and cross-stream migration 
of polymer away from the walls. Stress-induced polymer migration in a Couette 
device has been studied by Apostolakis et al. [1] using the two-fluid Hamiltonian 
model of Mavrantzas and Beris [IS] . The authors found that at large Deborah 
numbers polymer chains migrated from the outer towards the inner cylinder, 
in qualitative agreement with the theoretical predictions of Aubert and Tirrell 
[5] and the pioneering experimental results of Shaffer et al. [55] an d Dill and 
Zimm [12] showing DNA migration in high molecular weight DNA solutions 
between rotating concentric cylinders. In the present paper, the Kramers-type 
expression derived in Section 13.21 for the elastic stress takes into account the 
presence of walls and the consequent restrictions that this imposes on molecule 
configurations. Moving in a direction normal to the wall and into the fluid will 
generally mean that the configurations available to a dumbbell change and with 
this the evaluation of the stress tensor. It is the variation of stress in the normal 
direction to the wall that drives migration and leads to a non-uniform number 
density in the present case. An expression for the difference in the polymer 
and solvent velocities when the fluid surrounding a dumbbell is considered to be 
the solvent alone is given in Section [3.31 Although no explicit account is taken 
of hydrodynamic interaction in the developments of Sections [2] and |3j we show 
in Section 2] that considering the fluid surrounding a dumbbell to be the poly- 
mer solution leads to an expression for the perturbation to the solvent velocity 
that is proportional to the hydrodynamic force exerted by the dumbbells on the 
solution. 

Although the present paper is in the first place concerned with the derivation 
of a Fokker-Planck equation for dumbbells of FENE type in a non-homogeneous 
flow when bead inertia may be neglected, the underlying governing equations 
with the full account of the bead inertia are described in detail in Section [2] 
The derivation of the equations in the inertialess limit is performed in Section 
[3J In Section 13.11 we follow Degond and Liu [IT] in writing the full phase- 
space distribution function as a perturbation series about its inertialess value 
and show rigorously that the inertialess case is just the limit as the bead mass 
goes to zero of the full set of equations when inertia is included. The same 
result was obtained by Schieber and Ottinger J5B] who also concluded that 
both equilibration in momentum space and neglect of the acceleration term in 
the bead equations of motion are consequences of allowing the bead mass to 
tend to zero, and not independent assumptions as stated by Bird et al. [7J 
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and Curtiss et al. [5] . The two fluid theory of Doi [T3] , Doi and Onuki [T3] and 
Milner |19| has been influential in our approach to deriving the equation of linear 
momentum and of a Kramers expression for the elastic stress tensor in Section 
13.21 The solution is viewed as two coexisting continuous media, the polymer 
and the solvent, moving with different velocities. A separate momentum balance 
equations is written down for the solvent phase while the macroscopic average 
velocity of the polymer phase is obtained by averaging from the mesoscopic 
description. In our derivation we mostly suppose that the fluid surrounding 
a dumbbell in dilute polymer solutions should be taken to be the solvent and 
not the volume fraction average of the solvent and polymer, in contrast with 
some other kinetic theory descriptions of dilute polymer solutions [H [5] . We 
show, however, that assuming that the velocity of the surrounding fluid is the 
same as that of the solution as a whole is also possible in our framework, cf. 
Section |U The restriction of conformations available to model polymers near a 
solid wall mean that our derived Kramers-type expression for the elastic stress 
tensor implies, as mentioned above, a stress gradient in the normal direction to 
a wall and further, that off diagonal components (in particular, the elastic shear 
stress) are zero there. Sufficiently far from the wall, or in homogeneous flows, 
our Kramers-type expression simplifies to the usual one found in the literature 

tang. 

In reworking the kinetic theory of Bhave et al. [5] for non-homogeneous 
flow of dilute polymer solutions, Beris and Mavrantzas [3], referred to above, 
used truncated Taylor series about some spatial position vector r to obtain an 
expression for the polymer mass flux at r and then, using the polymer species 
conservation equation, wrote down an equation for the dumbbell number density 
which was common to all three formulations (5j UHl [23] considered. In Section l331 
we obtain the number density equation without resort to approximations and 
find that, if we chose the fluid velocity to be that of the solution rather than of 
the solvent, the exact number density equation would be precisely that derived 
in Beris and Mavrantzas's paper. This might explain why, despite differences 
in the constitutive equations appearing in [5]|18[f2^]. there is perfect agreement 
on the number density equation. 

In the homogeneous flow of a FENE-P or Oldroyd B fluid, substitution of 
the Kramers expression for the elastic stress tensor into the evolution equation 
for the dumbbell conformation tensor yields a closed-form constitutive equation. 
We seek to show that the kinetic theory presented here does not allow for such a 
closed-form constitutive equation in the case of non-homogeneous flows without 
resorting to closure or some other approximations, however. An alternative to 
a closure approximation would be to follow Bhave et al. [5] and simply employ 
the homogeneous form of the Kramers expression for the elastic stress tensor. 
This would, however, be inconsistent with retaining a diffusion term in the 
non-homogeneous Fokker-Planck equation, as we take pains to show in Section 

El 
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2 Kinetic theory 



We consider a dilute polymer solution filling a domain £1 C M. d . In the spirit 
of the two fluid theory of Doi [T3], Doi and Onuki [T3] and Milner [T^] we 
assume the solution to consist of two phases, polymer and solvent, with two 
distinct macroscopic velocity fields v p and v s , respectively. The description of 
the solution is completed by specifying the polymer volume fraction ip (and, 
consistently, the solvent volume fraction 1 — ip) which is allowed to vary in space 
and time. The solvent phase is assumed to be an incompressible Newtonian 
fluid and its velocity is governed by the macroscopic Navier-Stokes equations 
with the addition of the force term describing the interaction between the two 
phases. In contrast, the polymer phase is modelled on the mesoscopic level 
and the macroscopic quantities v p and ip are obtained through averaging. More 
specifically, the polymer molecules are represented by dumbbells consisting of 
two beads (more precisely, of two point masses), each of mass m, joined by a 
massless spring, as shown in Fig. [TJ The position vector of bead i (i = 1,2) is 
denoted by r,;(t) and its velocity by V = »V The velocity of the surrounding 
fluid, i.e. that of the solvent, at the point having position vector is denoted by 
Vi = v s (r.i). The equations of motion of the ith bead may therefore be written 
down as a first-order system in the following way: 



where £ denotes a drag coefficient, F\ is the spring force acting on bead i and 
Wi, W2 are two mutually independent d-dimensional Wiener processes. We 
assume that F\ = —F% = F, where F is a function of the dumbbell end-to- 
end vector T2 — T\ alone, ks and T are, in the usual notation, the Boltzmann 
constant and the (absolute) temperature. 

Remark 1 The velocity Vi = V s (ri) in (0J) should be understood as that of 
the solvent averaged over a box of size greater than that of a bead but smaller 
than the distance between the two ends of a dumbbell. We further assume that 
the fluctuations of the solvent velocity around v s can be disregarded over such 
a box. This means, in particular, that the mechanism of hydrodynamic interac- 
tion is not included in our picture, i.e. the fluctuations of the solvent velocity 
induced by the vibrations of a bead of a dumbbell are assumed to fade away on 
a length smaller then the distance to the other bead. We admit that this simpli- 
fying assumption of separation of scales is not necessarily satisfied in reality. A 
more adequate modelling may be achieved in the approach of Coupled Langevin 
Equations l HH. ]21^ where the solvent velocity is treated in essentially the same 
way as the polymer one, i.e. its random fluctuations on the length scale of a 
dumbbell are not neglected so that v s satisfies a stochastic (partial) differential 
equation. Such an approach is beyond the scope of the present paper and we 
content ourselves here with only a macroscopic description of the solvent, treat- 
ing only its averaged velocity. The choice of v s in the drag force in (QJ) is thus 
somewhat arbitrary. One could consider other options such as, for example, 



mdV t = (-((Vi - Vi ) + Fi) dt + y^Ck^TdW, 



(1) 
(2) 



dri = Vi dt, 



assuming that the averaged velocity of the surrounding fluid is that of the entire 
polymer-solvent solution. We reexamine this option in Section [^J 

We now introduce the 2d-dimensional vectors r = (rf, rJ) T , V = (V 1 T , V 2 T ) T , 
F = (F 1 T ,.F 2 T ) T , v = (vf,v%) T and denote by *(r, V,t) the dumbbell distri- 
bution function, defined such that ^(r, 1 ^,^ dr dV is the expected number 
of dumbbells having bead position vectors in the interval [r, r + dr] and bead 
velocities in the interval [V, V + dV] at time t. The Fokker-Planck equation 
equivalent to (d])-© may now be written down as 

^ + V r • (V*) + -V V ■ ((-C(V -v)+ F)*) = ^V 2 V *. ( 3 ) 

at m m 

Equation ([3]) should be equipped with some boundary conditions. We assume 
that the velocities V, can take any value in R d and <3/ vanishes sufficiently fast 
when Vi — > oo. We do not specify any particular boundary conditions for VP 
as r, approaches <9f2 for the moment. However, if there is no flux of the fluid 
across the boundary of f2, these boundary conditions should be chosen so that 

V 2 ■ n* da r2 dV = 0, Vri G fi, (4) 

an 

V x ■ n* da ri dV = 0, Vr 2 e 0, 

where d(T ri denotes an element of the surface as rj runs over 917, and the integral 
sign without a subscript (here and elsewhere) indicates the integration over the 
whole of R d or M. 2d . Indeed, take any volume u> inside fi and consider the 
dumbbells whose bead number 1 lies in uj. The beads numbered 2 in these 
dumbbells cannot leave the domain so that their flux through the boundary 
of fl should be 0, i.e. 

V 2 ■ n^> da r2 dridV = 0. 

Since this reasoning holds for any w C fi, the first equation in (j4j should be 
satisfied for all r\ € fi. The second equation in (Q| is established in the same 
way. 

We define the number density of dumbbells as 

N(x,t) := i (J jf *k=« dr 2 dV + J jf *k=«. dr.dV^j . (5) 

and the average polymeric velocity as 

v P : = ^ (/^ V^\ ri=x dr 2 dV + J J V 2 ^\r 2=x dndVj . (6) 

The polymer volume fraction tp is related to the number density N through 

<p = NV d , (7) 
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where Vd denotes the volume of a single dumbbell. 

It follows from the Fokker-Planck equation ([3]) that the number density TV 
should satisfy the continuity equation 

dN 

+ V x • (v p N) = 0. (8) 

To see this, let us integrate with respect to j*j over to get, for i = 1,2, 
respectively, 

— dr 2 dV + V ri • / / VL* dr 2 rfF = 0, (10) 
n or J J n 

where the divergence theorem in combination with |4]) has been used for sim- 
plification. Eqn. ([9]) is valid for any value of r 2 and Eqn. (flO|) for any value of 
r\ . We may therefore choose the value of r 2 in ([9]) and of J*i in (fTQ| to be equal 
to a;, add the two equations and divide by two. This gives us 



dN . 1, 
~W 



(J ^V 2 *| r2=a! dndV + J J^V^\ ri=x dr 2 dV^j =0, 
which is the same as ©. 



2.1 Continuity equations 

The continuity equations for the two fluid phases may be expressed as 



3(i - y) 

dt 



^+V (B -(^ P ) = 0, (11) 
+ V x • ((1 - tp)v a ) = 0. (12) 



Note that despite the similarity between (fTTj) and (|12p their mathematical role 
in our modelling is different: the polymer phase continuity equation pip fol- 
lows immediately from ([7|) and (JSj) while the solvent continuity equation (|12|) is 
independent and expresses the incompressibility of the solvent. Eqn. (|12| just 
means that the volume of solvent in any small material volume 5V(t), contain- 
ing at all times the same solvent molecules, remains unchanged even though the 
solvent fraction 1 — (p in the same material element may vary with time. The 
changes in if are due to the fact that as SV(t) is tracked in time, polymer is con- 
tinuously streaming through the surface of 6V(i) (because, in general, v p ^ v s ). 
It is in the sense that the volume occupied by the same solvent molecules is time 
invariant that the solvent (and, indeed, the polymer) is incompressible. 

The continuity equation for the volume-averaged velocity of the polymer- 
solvent solution (it, say) may be computed from the volume fraction average 

it = (pv p + (1 - cp)v s . (13) 
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Adding the two continuity equations (|TT|) - ([T^T) . we see that the two-phase fluid 
is incompressible: 

V x -u = 0. (14) 
2.2 Linear momentum equations 

We turn now to the elaboration of the linear momentum equations for the 
solvent and polymer phases. Considering a material volume SV(t) of a size 
greater than the size of a bead but smaller than the distance between the two 
ends of a dumbbell, we recognise 

d f 

- p s (l-<p)v s dV, (15) 

at JSV(t) 

as the rate of change of momentum of the solvent phase in SV(t), where p s 
denotes the (constant) solvent density. Equating (|T5|) with the sum of forces 
acting on the solvent in SV(t) and use of the Reynolds transport theorem and 
the solvent continuity equation (fT2j) yields, in the usual way (since SV(t) is 
arbitrarily chosen) that 

Dv 

p s (l - = Vs^lVs + Vs^x^x ■ v s ) - S7 x p s + f. (16) 

where the first two terms term on the right-hand side of (fTB"|) are due to viscous 
forces, p s is the pressure in the solvent and / = f(x, t) is the density of the force 
exerted at a point x and at time t on the solvent by the dumbbells. The material 
derivative here is associated with the solvent velocity, i.e. = + v s ■ V x . 
The key observation in defining / is that all the interaction forces between the 
polymers and the solvent are contained in the equation of motion of dumbbells 
(fTj). Leaving to one side the spring force there (which acts inside the polymer 
phase), the action of the solvent on the polymer is described by the friction 
force — C(Vi — Vi) (on the bead number i) and the Brownian force. We argue 
furthermore that the Brownian forces do not contribute to /. Indeed, they are a 
part of the mesoscopic description of the fluid and they account for random kicks 
experienced by the polymer molecules from the surrounding solvent molecules. 
The quantity / on the other hand, exists on the macroscopic level, hence its 
computation involves averaging on volumes that are sufficiently large compared 
to the micro-scale. The random kicks should cancel each other during this 
averaging since the probability of any such kick is the same as that of one in the 
opposite direction. This observation is also consistent with the mathematical 
fact that the expectation of an increment of the Wiener process during any 
time interval is 0. Therefore, only the friction (hydrodynamic) forces should be 
incorporated into the formula for /. By the principle of reciprocal action we 
thus obtain the following expression: 

/(*,*)= / J C(Vl-«i)*k=« dr 2 dV + J J ((V 2 -v 2 )y\ r2=x d ri dV, 

(17) 
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where the first integral accounts for the dumbbells whose bead number 1 is at 
position x and the second integral accounts for those dumbbells whose bead 
number 2 has position vector x. By comparing (fTT)) with the definitions and 
^ of the number density and the polymeric velocity and noting that Vi = v s (x) 
when Ti — x we easily obtain another useful formula for the force /: 

f = 2N((v p -v s ). (18) 

The description of our model of the solution of dumbbells that takes into 
account the inertia of the beads is now complete. The system of governing 
equations thus consists of the Fokker-Planck equation §3§ for the polymer phase 
and the Navier-Stokes-type equations (|T2"j) - (jT(3)) for the solvent phase coupled 
through the expressions for the number density, volume fraction and inter-phase 
force density, i.e. (JSJ), ([7]) and (fTT)) . respectively. We emphasize that these 
equations do not involve the averaged polymer velocity v p which, together with 
the total velocity it, should be regarded in our framework as a quantity emerging 
in a post-processing of the basic equations. Indeed, v p can be always computed 
via once the distribution function *B is known. In particular, we don't 
really need a linear momentum equation for the polymer phase. However, it 
is interesting to write down such an equation in order to compare it with the 
corresponding equation for the solvent phase (fT6| and to see the relations of 
our modelling to the previous work of other authors. Assuming that the fluid 
fills the whole space or that *B vanishes on dSl, this equation reads (see the 
derivation in Appendix B) 

f dv \ 

Pp(p [ ~~dt + Vp ' WxVp ) = ~ 2N ^ v p ~ Vs "> + Vx ' f ' ~ Pp Vx ' (War(V)), (19) 

where p p := 2m/Vd is the polymer density, and Var(V) is the variance of V 
defined by 

Var(V) = ^(7 [ (Vi - v p )(Vi - v p )*\ ri = x dr 2 dV 



n 



in 

The tensor -P in (11911 is defined as 

.1/2 
-1/2 



(V r 2 -»p)(V 2 -tv)*U=« d '" 1 dV'). (20) 



1 /2 

(x,t)= [ [ [ H>\ ri = m+u _ i)q qF(q) dVdqds (21) 

J -1/2 J J , .,,,,!,„ 



r 2=a! + (s+i)q 



This is a modification of the Kramers expression for the polymeric contribution 
to the stress tensor, in which only the inter-bead spring force is taken into ac- 
count. A similar tensor arises t s will arise in our study of the vanishing bead 
inertia limit, cf. (|Tf|) and Appendix A. The conceptual difference between the 
linear momentum equation for the solvent (|16|) and that for the polymer (|17p 
resides in the presence of Var(V) in the latter, which measures the fluctuations 
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of the polymer velocity, while the fluctuations of the solvent velocity are disre- 
garded. This is consistent with our basic assumption that the solvent can be 
wholly modelled on the macroscopic level while keeping a mesoscopic descrip- 
tion for the polymer, as explained in Remark 1 on p. [SJ Momentum balance 
equations for the two phases in the two fluid models of Doi and Onuki [T3] and 
of Milner [H] , analogous to those in (|TH)l and (TT5]) , may be found in Equations 
(4.12)-(4.13) and (14) of these papers, respectively. Note that their description 
of both phases is macroscopic which manifests itself in the absence of the terms 
like Var(V). 



3 Negligible bead inertia 
3.1 The Fokker- Planck equation 

We now proceed to obtain the Fokker-Planck equation for the dumbbell distri- 
bution function satisfied in the limit of negligible inertial forces. To this end, 
and following Degond and Liu jll] , we use the following scaling and change of 
variables: 

m = £ 2 , p = eV. (22) 

To motivate this scaling, we note that the energy equipartition theorem yields 
(mV^ = 3ksT in thermal equilibrium so that, supposing that the velocity 
distribution is not far from the equilibrium one, we get that the characteristic 
value of the velocity should be of order yJh^Tjm, so that V behaves like 1/e 
as e — > 0. 

Writing out a regular perturbation series for \& 

* = *o + e*i+e 2 *2 + ..., (23) 

we are looing for a closed-form equation for V^o in which the contributions of 
order e and higher are neglected^ After the change of variables (|22p. the Fokker- 
Planck equation ([3]) becomes 

9* 1 1 1 

— + -V r ■ (p*) + -V p ■ ((C« + F)V) = ? CQW, (24) 

where 

Q(tf ) := V p ■ (ptf ) + fc B TV p *. (25) 
Substituting the perturbation series for ([23| into (|24| leads to 

+ V r ■ (ptt fc ) + V p ■ (((v + F)* fc ) = CQ(*fc+i), (26) 

1 Treating e as a small parameter is strictly speaking meaningless since it is a dimensional 
quantity. However, by performing the standard and well-established non-dimensionalization 
of all the variables we can easily see that the non-dimensional counterpart of e is \J Xb / Ah 
where Xb = m/C ls the characteristic time scale for the velocity fluctuations and A^f = C/4.ff 
is the characteristic relaxation time of the dumbbell, H being the scale for the spring force. 
The condition Xb << Ajy is fairly well satisfied in most experiments, cf. 1251 1261 . Here, we 
keep the dimensional parameter e only to make the notation more readable. 
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for k = — 1, 0,1,.. ., where we understand and ^ -2 to be zero. We introduce 
the marginal distribution functions ipk{v,t) as 

Mr,t)= J * fe dp, (k = 0,1,2,...). (27) 

Then, integrating (|26p throughout with respect to p over all p — space we arrive 
at 

^ = -Vr • / P^k+i dp, (k = 0,1,2,.. .). (28) 

We first set k = -1 in which yields <5(^ ) = 0. Now, from (g5]) it is 
not hard to see that Q may be rewritten as 

Q(f ) = fc B TV p • (MV P ) , (29) 

where M is the Maxwellian function 

M(p) = Ccxp(-p 2 /2k B T), (30) 

C is a normalizing constant (chosen to make the integral of M equal to 1) 
and p 2 :— WpW 2 - From the divergence form of operator Q in (|29l) . it may be 
shown that the kernel of Q is the linear span of the Maxwellian function M(p) 
defined in (|30|) . cf. [11]. Consequently, ^o(r,p,t) for an y fixed r and t should 
be proportional to M(p). In view of ([2T)l and the normalization of M(p), this 
means that ^(^P^) = ipo( r ,t)M(p). 

The remainder of the task of obtaining a Fokker-Planck equation for ipo will 
be occupied with deriving, when k = 0, an expression for the right-hand side of 
(|28p in terms of "!/V Note, therefore, that use of the identity 

p+-^-VpM = 0, (31) 
and repeated application of integration by parts gives for any \f 

J pQ(tt) dp = k B T J pV p ■ (mV p dp = -fc B T y MV P f |Q dp, 

= fc B T / ^V P M dp=- J pt> dp. (32) 

Eqn. (|26p provides another way of deriving an expression for the first integral 
in (|3"2"|) when ^ is replaced by ^k+i- multiplying (|2l)|) throughout by p/£ and 
integrating with respect to p yields 

pQ(^ fc+ i) dp = \ ( ^- [ p* fc _i dp + V r ■ [ pp^ k dp 



I 



pV p • ((C« + F)tf *) dp 
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Performing an integration by parts on the last term here and combining it with 
(|32|) we arrive at 

- J pffc+i dp = i J p^k-i dp + V r • J pp^ k dp - {(v + F)tp k ^ . 

(33) 

Now set k = in (|33|) . As remarked already. = 0. Moreover, using the 
identity (|31| and integration by parts we see that 



J pp^o dp = Vo J PPM dp = -k B T^ J p (V P M) dp 

= k B T^ 8 J M dp = k B Tip S. (34) 

Therefore, (f33|) with k = is simplified to 

J p^ 1 dp = 1 ((C« + J> - fc B TV r Vo) • (35) 

From this, setting fc = in (|28p we get the closed equation for ipo which is now 
a function of only r and t. This is the well-known Fokker-Planck equation for 
a dilute solution of dumbbells in the inertialess (e — > 0) case: 

^ = Jv r ■ (fc B TV^o - (C« + Ftyo) ■ (36) 
The right-hand side of (|3"6")) is just a compact notation for 
1 2 

- ^ V ri ■ (k B TV r ^o - (C«i + ■ 

The Fokker-Planck equation may be recast into a more familiar form in terms 
of the position vector of the centre of mass, x, and the end-to-end vector, q, 
defined as 

x := -(n +r 2 ), q := r 2 - ri. (37) 

For clarity of notation we will denote the function ?/>o of n, r 2 and i by ip 
when it is considered as a function of x, q and i. We also recall the convention 
F\ = —F% = F = F(q). Performing the change of variables from r\, r 2 to x, q 
we thus arrive at 

dtp _ / 2k B T . 2F . , 

= V q • ( ~ £~ v ^ + — V + («! - w 2 )V 

^V^-f^]^]. (38) 
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3.2 Polymer force acting on the solvent and the elastic 
stress calculator 



We now turn to the calculation of the force vector / appearing in ([16)) in the 
case of negligible bead inertia. By writing Vi = Pi/y/m (= Pi/e) and retaining 
only the O(e ) terms in (fl~7|) . the force vector is given by 

f(x,t) = ( / / Pi*i|n=a! dr 2 dp - ( / vii/) \ ri = x dr 2 
J Jo Jo 

dr idp -( v 2 ipo\r 2 =x dn. (39) 
J Jn Jn 

The first moments of can be related to ipo as shown in equation (|35|) . Substi- 
tuting (|35p into ([321) I we see that / can be rewritten as a sum of two components 
f s (due to the spring) and /* (due to thermal motion) 

f(x,t) = f s (x,t) + f(x,t), (40) 

with 

f s (x,t) = / FiV>o|n=* dr 2 + / F 2 ^ \ r2=x dr u (41) 
Jn Jn 



f{x,t) = -k B T I I V ri ^o|fi=a! dr 2 + / V T . 2 V'o|r 2 =x dri 
\Jn Jn 



(42) 



In Appendix A it is shown that the component f s can be represented as the 
divergence of a tensor / s = V x ■ t s with t s given by 

T s (x,t)= f f ij)(x + sq,q,t)qF dqds, (43) 

J -1/2 J Q s (x) 

where Q s {x) is the set of vectors q such that x + (s ± l/2)q G f2. For the other 
component /' of the force, we immediately obtain 

f(x,t) = -k B TV x / Vol T*i — X dr 2 + / Vol dn = -2k B TV x N, (44) 
\Jn Jn J 

where we have used the formula for the dumbbell number density in the iner- 
tialess case 

N(x,t) = \(J Volr^x dr 2 + J Vo|r 2 =x dr^j , (45) 

which is a direct consequence of the general definition Putting (|4"3"|) and 
together we conclude that 



r(x,t) := [ ( ip{x + sq,q,t)qF dqds - N(x,t)k B TS (47) 

J -1/2 JQsM 



f = V x • t s - 2k B TV x N = V, • t - V xPp , (46) 
where, we have introduced the elastic stress tensor as 

.1/2 
'-1/2 JQ a (x) 
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and the polymeric contribution to the pressure as 

p p (x,t) := N(x,t)k B T. (48) 

The reason for decomposing t s into a sum of r and a pressure-like term is to 
ensure that the elastic stress t vanishes at equilibrium, consistent with classical 
definitions, as in [7J, for example. 

Taking Eqns. (|16p and (|46p together we may now write down the equation 
of linear momentum for the solvent phase as 

D v 

Ps {l - <p)-j£ = Vs^lvs + VsVx(Vx ■ v a ) - S7 xP + V x • r, (49) 
where p := p s + p p is the total pressure. 



Some remarks on (|47|) 

1. If we neglect the presence of the walls by writing Q s (x) = Q, independent 
of x and s, for any x and s, the subtraction of the isotropic term on the 
right-hand side of the expression (|47p for the elastic stress tensor ensures 
that in equilibrium (i.e. when \7v = 0) t = 0. See, for example Eqn. 
(13.2-18) of [7J. 

2. In the non-homogeneous case, where full account is taken of wall effects, 
we note from (|47p that the elastic stress tensor r will not, in general, 
vanish in equilibrium. Moreover, on smooth boundaries dfl which are 
impenetrable by the dumbbell beads the configuration space Q s (x) has 
zero measure in d dimensions, leading to the boundary evaluation of the 
elastic stress 

r(x,t) = -N(x,t)k B TS Va; g dQ. (50) 



3. An expression similar to (|47|) for the elastic stress was derived by Biller 
and Petruccione in H>1 [5U and used in [TB] . 

4. We can provide an alternative approximate formula for the elastic stress 
that does not involve the auxiliary variable s. To this end, we recall 
the second line in (|69p . which is a weak expression for f s valid for any 
compactly supported test function g, and transform it by expanding g(x± 
q/2) in a Taylor series about x and integrating by parts: 

f s gdx = -[ [ Ftp- (q-V x g + ^-(q-S7 x ) 3 g + ■■■) dqdx 

n JnJQ(x) V ^ 4 / 



/ g -v x - 


\[ qFi'dq) 


In 


\JQ( X ) J 



1 

24 



g V x \ / qF{q-V x fi>dq dx 
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We thus see again that the vector f s is represented as a divergence of a 
tensor field. Invoking the same expression for f t as before and subtracting 
the isotropic term we obtain that the elastic stress introduced in (|4"T)) can 
be also approximately computed via 

t(x, t)= I qFi/j(x, q,t)dq+ — [ qF (q ■ V x ) 2 ip(x, q,t) dq + ■ ■ ■ 

Jq{x) ^ 4 Jaw 

-N(x,t)k B TS. (51) 

Retention of only the first term in the Taylor series in (|51[) would then 
yield the usual Kramers expression for the elastic stress tensor: 

r = ( qF ) - Nk B TS, (52) 

where we have used the notation ( • ) for the ensemble average 

(•):= f -TP(x,q,t)dq. (53) 
JQ 

3.3 The polymer velocity and the number density equa- 
tion 

The polymer phase velocity, v p , defined by (|6]), can be easily computed in the 
inertialess limit thanks to the formula (|18p . Indeed, substituting the expression 
for / from we arrive at 

1 :f = v s + -5-V. • r - ^V X N. (54) 



2iVC 2(N 2(N 

Inserting this expression for v p into (|8|) we finally arrive at the equation for 
the polymer number density in the inertialess limit 

f +V,- (v s N) = ^VlN ±V X V X : r, (55) 



where we use the notation V X V X '■ t in (|55|) to mean the divergence of the 
divergence of r, i.e. 

i=i j—i j 

Note that, unlike Beris and Mavrantzas [J], we have not needed to use Tay- 
lor's theorem and the neglect of higher-order terms in order to arrive at the 
convection-diffusion equation ([S3)) for N. The exactness of this equation might 
explain, therefore, why it was that the consistent derivation of the non- homogeneous 
kinetic theory, the continuum two-fluid Hamiltonian model and the body-tensor 
continuum formalism considered by these authors resulted in identical descrip- 
tions of the equation satisfied by the polymer number density. 
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3.4 Summary of the model with neglected bead inertia 

As the governing equations of the model in the inertialess case are scattered 
throughout the article, it can be helpful to recapitalize them here as follows: 
the solvent velocity v s is obtained from the equations of motion (|4"9")l and (fT2"j) , 
where the elastic stress tensor r in (j4"9")l is found from Eqns. (|4"7} -(j4"5" j) and ip 
in (j4"Tj) is the solution to the Fokker-Planck equation with V{ = v s (r.i), 
while the polymer volume fraction ip is calculated from ip through (|45|) and J7J) . 
The continuity equation (|12[) can be reaplaced by an equivalent equation for the 
divergence of v s : 

Iv. Vs = Mv2jV-iw:T, 
V d 2C 2C 

which is easily derived with the aid of (f5"5"j) . If needed, the average polymer phase 
velocity v p can be computed then by (|54[) and the volume-averaged velocity of 
the solution it by (fT3|) . 



4 A modification of the theory with the assump- 
tion "surrounding fluid = solution" 

All of the preceding theory has been developed assuming "surrounding fluid = 
solvent" in the drag force term of ((T]). As has already been mentioned in Remark 
1 on p. [S] this hypothesis is rather arbitrary although it seems appropriate in 
very dilute solutions and in the absence of hydrodynamic interaction. It is also 
the underlying assumption in the works of |13l 1141 119] . However, we may also 
consider an alternative hypothesis "surrounding fluid = solution" in ([1]), which 
can be seen as a very crude way to account for the hydrodynamic interaction. 

The theory of Sections O [3] may be easily modified to accommodate this 
case: in the dumbbell equations of motion ([I])-©, for example, the only change 
is to set Vi = u(r.i), i = 1,2, u being the averaged solution velocity. The kinetic 
theory for the polymer phase remains essentially the same and the Fokker-Planck 
equation (|38p still holds in the inertialess limit, again assuming Vi = u(ri). 
However, some changes should be made to the coupling of these mesoscopic 
equations and the macroscopic motion equation of motion for the solvent (|16[) . 
In the equations below, and for ease of notation, let us keep the definition of 
the force / to be the friction force 2N((v p — v s ) between the polymer and the 
solvent, as it was in the preceding section. 

Without entering into all the details of the derivation, which are essentially 
the same as in the preceding section, we just summarize here the governing 
equations stemming from the hypothesis "surrounding fluid = solution". They 
include the Fokker-Planck equation (f3"5|) with the substitution V\ = u(x—q/2, t), 
v-2 = u(x + q/2, t), while the equations for the macroscopic quantities v s , u and 
r become 

Dv 

p s (l - <P)~^ = Vs^lvs + VsVx(Vx ■ v s ) - S7 x p s + f, (57) 



1G 



) + V* • ((1 - <p)v s ) = 0, 



dt 

.1/2 



r(x,t)= [ [ ip(x + sq,q,t)qF dqds - 2N(x,t)k B Td, 

J -1/2 JQ s (x) 

u = Vs + K^V) Vx ' T ' (58) 

with iV calculated from (being careful to replace v 8 in this equation with 
u). Since {v p — u) = (1 — <p)(v p —v s ), the force / on the right-hand side of (|57|) 
can now be calculated from 

1 - </5 

Note, therefore, that the force between the polymeric and solution phases in 
Eqn. (|57p is no longer the divergence of a tensor because of the division by 
(1 — tp). For the same reason, we are not easily able to interpret a part of the 
elastic stress as a pressure-like term. 

By writing the solution velocity u as the sum of the solvent velocity v s and 
a disturbance velocity v' s , we see from and (|58|) that 

< = Y C f = W^) c(Vp " u) = W^) c(Vp " {Vs + <}) - (59) 

Although Eqn. (|59| holds as a direct consequence of the definition (fT3|) of 
the solution velocity u whether we consider the surrounding fluid to be the 
solvent or the solution, the interpretation is rather different in the two cases. By 
"surrounding fluid" we really mean the fluid in the immediate neighbourhood of 
a dumbbell. In Section EOl v' a would just be the difference between the solution 
velocity and the solvent velocity, although the velocity field in the immediate 
neighbourhood of any individual dumbbell is considered to be affected by other 
dumbbells only in so far as they contribute to the elastic stress in the equation 
of motion (|49|) . When the surrounding fluid is considered to be the solution, v' s 
is the actual disturbance to the velocity field in the immediate neighbourhood 
of a bead caused by the other dumbbells beyond their contribution to the elastic 
stress in the solvent equation of motion. From (|55)) , v' s is seen to be the product 
of a function and the force exerted by the polymer phase on the solution and is 
much simpler than in the case of full hydrodynamic interaction. 

5 Non-dimensionalization of the governing equa- 
tions and elaboration of important limiting cases 

Let V and L denote characteristic velocity and macroscopic length scales, re- 
spectively, and N av be a (space- and time-) averaged value of N. The homo- 
geneous version of the Fokker-Planck equation (|38|) for a Hookean dumbbell 
in equilibrium (v s = 0) tells us that a characteristic mesoscopic length scale 
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la := \J tx(qq) jd satisfies i^H = k B T where H is the spring force constant for 
a Hookean dumbbell. 

We may introduce non-dimensional variables: 

* = i'"- = = t>* = ^' 9 =f' F (60) 

This rescaling leads to the non-dimensional Fokker-Planck equation, cf. (|38[) . 

^! = V„,- f _^v„-v* + — ^*-|(9*-V IB .)«! + — f-^ (9*-V„.) 3 « 
Di* 9 \2De q v 2De V \ KH xjs 2A\LJ KV ' 

+ v,. • ( -L. V v..*- ~ { 7 ( I V (<** " V- ) 2 < + • • ■ L* V 



8L>e U/ 4 V L 



(61) 



where 



De:=§- L , (62) 

is a Deborah number and the non-dimensional distribution function ip* is intro- 
duced so that ip = ip*N av /£Q. The terms discarded in (|6"Tj) are of order 4 and 
higher in £q/L. 

In order to give a reasonable scaling for the elastic stress we remark that it 
is of order N av kBT(^/ (4H)(Wv + Wv T ) in the regime of small deviations from 
equilibrium (in fact, this is the first term in the Taylor series for r as \7v — > 
in the case of Hookean dumbbells). The dimensionless elastic stress tensor r* 
is thus introduced as 

(63) 



T 



N av k B TDe 
and the expression (fSTj) for r becomes 



r*(x*,t*) = —\ [ q*F*i)*(x*,q*,t) dq* 

n 

[ q*F*(q*-V 3! .) 2 i>*(x*,q*,t)dq* + 



1 

24 
N*S 



(64) 



where we have written out only those terms up to O(£o/L) 2 , as in (j6~T]) . The 
expression (j54")l for the polymer phase velocity becomes 

'~V fv«..T*--^V..i\rY (65) 



v s 8N* \ L ) V De 
The dimensionless number density equation (|55[) now reads 

^ + V,. • KAH = I (|) 2 (-L-VI.N- - V..V,. : r« ] . (66) 
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Finally, the characteristic total viscosity may be introduced as 77 = rj s +N av kBT(/(4:H) 
and the corresponding Reynolds number is Re = p s VL/r). This leads to the 
non-dimensional equation of motion 



Re(l - <p) 



Dv* 



(67) 



Dt* 



where 77* and 77* are dimensionless solvent and polymer viscosities defined, re- 
spectively, as 77* = r)*/r) and 77* = N av k B T(/(4:Hr]). 

5.1 io/L negligible 

Let us first consider the equation set to be solved in the case of £q/L negligi- 
ble (known as locally homogeneous flow). In this case, there is no distinction 
between the phase velocities v* = v* = u* and the Fokker-Planck equation 
(|61|) assumes its standard form, in particular there is no diffusion in x there. 
Moreover, as remarked already in Section [XU and as may be seen from (|M|) . the 
Kramers expression holds for t* 



where the average is taken over the configuration space Q* , independent of 
x* . As is well known, a closed- form differential equation for (q*F*)* may 
be obtained in the case of FENE-P or Hookean dumbbells, for example, by 
multiplying the Fokker-Planck equation ([6"Tj) throughout by q*F* and using 
integration by parts. This equation, together with the equations of motion (|67jl 
for it* = v*, the incompressibility equation V • u* = 0, and the number density 
equation DN* /Dt* = 0, now form a solvable system, to be solved subject to 
suitable initial and inflow conditions. 

5.2 Non-homogeneous flows 

The major impediment to elaborating a system of equations to be solved in 
the non-homogeneous case, even when the dumbbell spring is Hookean, is that 
retention of terms of 0(£q/L 2 ) in (|B"Tj) , (fM| and (|6"B"|) means that the system 
cannot be closed even if one persuades oneself that the 0(£q/L 2 ) terms in the 
Fokker-Planck equation (|6ip may be neglected when the length scale for velocity 
variations is comparable to £q. See the arguments used by Beris and Mavrantzas 
|4j, for example. This is because, in contradiction to the form of the Kramers 
expression assumed by Bhave et al. [S] and Beris and Mavrantzas [3], for 
example, one must, for the sake of consistency with what is done elsewhere, 
retain all terms up to 0(£q/L 2 ) in the expression for t. On the other hand, 
if we are prepared to sacrifice consistency and to keep the terms of order (^0 / L) 2 
only in the x-diffusion term in ([61]) together with (|65|) and (|66|) . a closed- form 
constitutive equation for r can be derived in our framework as well. 



(q*F*)* - N*S, 



(68) 
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Conclusions 



We have presented a framework to model the dilute polymer solution in the spirit 
of the two fluid theory. Unfortunately, the obtained equations do not reduce 
to simpler constitutive equations in the case of strongly non-homogeneous flows 
in any self-consistent way, even if Hookean dumbbells are used to represent the 
polymer molecules. One should hope, however, that their numerical simulation 
is possible by stochastic or Fokker-Planck-based methods as is demonstrated, for 
example, in [16] for a similar model. We believe that the use of a perturbation 
series in going from the equations of motion with bead inertia to the inertialess 
limit provides a pleasing and rigorous method of deriving both the Fokker- 
Planck equation and a generalized Kramers-type expression for the elastic stress 
for this case. Moreover, this approach may be pursued further in keeping more 
terms in the perturbation series. This will provide not only additional terms in 
the momentum balance equation but also in that of the energy equation. 



A The averaged spring force in the inertialess 
case 

Multiplying the formula (pTTj) for the averaged spring force f s by an arbitrary 
vector- valued test function g vanishing on the boundary d£l, integrating over Q 
and changing the variables in the integral from dr\dr 2 to dxdq, we see that 

/ f s {x,t) ■ g(x) dx = / / FiV> -g(ri) dr 1 dr 2 + \ \ F 2 ip ■ g(r 2 ) dr 1 dr 2 
Jn Jn Jn Jn Jn 

= f [ Fi> ■ (g(x-q/2)- g(x + q/2)) dqdx (69) 
Jn J q{x) 

= — I / Fip(x,q,t)- / q -Vg(x — sq) dsdqdx, 

JnJQ(x) J -1/2 

where Q(x) is the set of vectors q such that x±q/2 <G fl. By making the change 
of variable x — > x + sq, we arrive at 



[ f s (x,t) ■ g(x) dx = - [ [ [ tjj(x + sq,q,t)F -(q-Vg)(x) dqdsdx, 
Jn Jn J -i/2 Jq 3 (x) 



where Q s {x) is the set of vectors q such that x + (s ± l/2)q € il. Using 
integration by parts, this leads to 

f s (x,t) = V x ■ f 1 [ iJj{x + sq,q,t)qF dqds. (70) 

J -1/2 JQ B (x) 

We recognize that the right-hand side here is the divergence of tensor r s as 
defined in (|gg]>. 
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B Derivation of the averaged momentum equa- 
tion for the polymer phase ( 1191 ) 

Let us first note that the definitions of the number density (JSJ), the average 
polymeric velocity © and its variance (12"0)) can be simplified thanks to the 
symmetry *(n, r 2 , V\, V 2 , t) = *(r 2 , n, V 2 , Vl, t): 

N(x,t) = J Jj\ ri = x dr 2 dV, v p (x,t) = -^ J J Vi*| ri=sc dr 2 dV, 

Var(F) = 1 y y (Vx - u p )(Vi - Wp)*^^ dr 2 dV. 
This leads immediately to the following relation 

Var(F) = 1 y y Vi*!* | ri=a! rfr 2 rfy - « p « p . (71) 

Let us now proceed to the derivation of (|19p starting from the Fokker-Planck 
equation ([3]). We multiply (J3J) by V\ and integrate throughout with respect to 
r 2 and V . In doing this we denote r\ as a; to be consistent with the formulas 
for N and v p above. After an integration by parts with respect to V, we thus 
obtain 

J J V ^ 1=x dr 2 dV + V X - J J ' ViVx^^dndV + J J ViVa-V,*.*^ 
= - / f(-C(V 1 -v 1 ) + F 1 )^\ ri=a: dr 2 dV. 



The third term on the left-hand side can be reduced to an integral over dfl 
with respect to r 2 and thus it vanishes if ^ is zero on the boundary dfl or, 
in particular, when is the whole space M. d . All the other terms in the last 
equation are simplified with the aid of definitions of N and v p and relation (|7ip : 

d NC 1 - 

g- t (Nv p ) + V x • {Nv p v p + NVm(V)) = — - v s ) + —f, (72) 

where f s is the averaged spring force 

f s (x,t) = J J Fi*| ri=a , dr 2 dV + y J F 2 V\ r2=x d ri dV 

= 2 [ [ F^\ ri=x dr 2 dV. 



The last quantity is almost the same as f s introduced by formula (|41[) in the 
study of the inertialess case. The only difference is that ipo featuring in the 
definition of f s is now replaced by ip = J ^>dV. We can thus repeat the 
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Figure 1: Elastic dumbbell consisting of two point masses of mass m, joined by 
a massless spring, x denotes the position vector of the centre of mass and q the 
end-to-end vector. 



reasoning of the preceding appendix B in order to rewrite f e as the divergence 
of a tensor 

f s = Vx • / / i>{x + sq, q, t)qF dqds := V x -f s , 

J -1/2 J 

where ?p is ipg rewritten in terms of r c and q. The definition of t s is rewritten 
in (|2ip in terms of "J. Multiplying (|72| by 2m, introducing the polymer density 
p p = 2m/Vd and recalling that N = (p/Vd, we now obtain 

d 

-7^(PpipVp)+\7 a .-(p p ipVpV p ) = ~2N((Vp-v s ) + \7 a .-T s -\7 x -(pp(pVa,r(V)). (73) 
Combining (|73|) with the continuity equation (jHJ leads to (fT9l) . 
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